#!/usr/bin/python
# -*- coding: UTF-8 -*-
"""
@author:Venus
@file:2100832-唐子钊.py
@time:2021/12/13/16：54
"""
import sys

sys.path.append(r'c:\programdata\anaconda3\lib\site-packages')
import numpy as np
import struct
import prettytable as pt
import time
import tqdm

# 训练集文件
train_images_idx3_ubyte_file = r'C:\Users\Venus\Desktop\Pycharm_Workspace\venv\智能分析与决策\数据集\mnistdataset\train-images.idx3-ubyte'
# 训练集标签文件
train_labels_idx1_ubyte_file = r'C:\Users\Venus\Desktop\Pycharm_Workspace\venv\智能分析与决策\数据集\mnistdataset\train-labels.idx1-ubyte'
# 测试集文件
test_images_idx3_ubyte_file = r'C:\Users\Venus\Desktop\Pycharm_Workspace\venv\智能分析与决策\数据集\mnistdataset\t10k-images.idx3-ubyte'
# 测试集标签文件
test_labels_idx1_ubyte_file = r'C:\Users\Venus\Desktop\Pycharm_Workspace\venv\智能分析与决策\数据集\mnistdataset\t10k-labels.idx1-ubyte'


# 解析图片数据
def decode_images(idx3_ubyte):
    # 读取二进制数据
    bin_data = open(idx3_ubyte, 'rb').read()
    offset = 0
    fmt_header = '>iiii'
    # 魔数          图片数       图片行数   图片列数
    magic_number, num_images, num_rows, num_cols = struct.unpack_from(fmt_header, bin_data, offset)
    image_size = num_rows * num_cols
    offset += struct.calcsize(fmt_header)
    fmt_image = '>' + str(
        image_size) + 'B'
    images = np.empty((num_images, num_rows, num_cols))
    for i in range(num_images):
        images[i] = np.array(struct.unpack_from(fmt_image, bin_data, offset)).reshape((num_rows, num_cols))
        offset += struct.calcsize(fmt_image)
    return images


# 解析标签数据
def decode_labels(idx1_ubyte):
    bin_data = open(idx1_ubyte, 'rb').read()
    offset = 0
    fmt_header = '>ii'
    magic_number, num_images = struct.unpack_from(fmt_header, bin_data, offset)
    # print("1:", type(num_images))
    offset += struct.calcsize(fmt_header)
    fmt_image = '>B'
    labels = np.empty(num_images)
    for i in range(num_images):
        labels[i] = struct.unpack_from(fmt_image, bin_data, offset)[0]
        offset += struct.calcsize(fmt_image)
    return labels


def show(labels, result, name):
    TP = FN = FP = TN = 0
    for i in range(len(labels)):
        if labels[i] == 1 and result[i] == 1:
            TP = TP + 1
        elif labels[i] == 1 and result[i] == -1:
            FN = FN + 1
        elif labels[i] == -1 and result[i] == 1:
            FP = FP + 1
        elif labels[i] == -1 and result[i] == -1:
            TN = TN + 1

    print(name)
    print("          y_hat = +1            y_hat = -1                  Total")
    print("y=+1     ", "TP=", TP, "              FP=", FP, "               TP + FN = ", TP + FN)
    print("y=-1     ", "FN=", FN, "              TN=", TN, "             FP + TN = ", FP + TN)
    print("All      ", "TP + FP =", TP + FP, "        FN + TN =", FN + TN, "       TP + FP + FN + TN = ", TP + FP + FN + TN)
    print('错误率=', ((FN + FP) / (TP + FP + FN + TN)) * 100, '%\n')

if __name__ == '__main__':
    start_time = time.time()
    # 获取数据
    print("-----1.载入数据-----")
    train_images = decode_images(train_images_idx3_ubyte_file)
    train_labels = decode_labels(train_labels_idx1_ubyte_file)
    # [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
    test_images = decode_images(test_images_idx3_ubyte_file)
    test_labels = decode_labels(test_labels_idx1_ubyte_file)
    # 数据预处理
    print("-----2.数据预处理-----")
    # 1.找到493维像素位置
    # 训练集
    train_images_pre = train_images.reshape(60000, 784)
    train_labels = train_labels.reshape(60000, 1)
    # 测试集
    test_images_pre = test_images.reshape(10000, 784)
    test_labels = test_labels.reshape(10000, 1)
    # 如果图像为0 则取1，否则取-1
    for i in range(len(train_labels)):
        if train_labels[i] == 0:
            train_labels[i] = 1
        else:
            train_labels[i] = -1

    for i in range(len(test_labels)):
        if test_labels[i] == 0:
            test_labels[i] = 1
        else:
            test_labels[i] = -1

    # 用0范数统计至少在600张图片中不为0的像素，作为特征提取
    # 一共有28*28=784个特征点，提取特征后发现有用特征为434个
    # 统计出符合要求的特征点数量
    feature = np.zeros([784, 1])
    for i in range(784):
        flag = np.linalg.norm(train_images_pre.T[i, :], ord=0)
        if flag >= 600:
            feature[i] = 1
    row = int(feature.sum())
    # 训练集
    train_images_feature = np.zeros([row, 60000])
    flag2 = 0
    index = []
    for i in range(784):
        if np.linalg.norm(train_images_pre.T[i, :], ord=0) >= 600:
            train_images_feature[flag2, :] = train_images_pre.T[i, :]
            flag2 += 1
            index.append(i)
    # 测试集
    test_images_feature = test_images_pre.T[index, :]
    # 首先要在开头插入全为1 的列，组成494 维行向量，在水平方向上平铺
    A = np.hstack([np.ones([60000, 1]), train_images_feature.T])  # 训练集
    b = train_labels
    # 2. 步骤1：QR分解
    print('-----3.进行QR分解-----')
    q, r = np.linalg.qr(A)
    print('-----4.完成QR分解-----')
    print('-----5.最小二乘计算-----')
    # 3. 步骤2 计算QTb
    QTb = np.dot(q.T, b)
    # 4. 步骤3：求解x
    x = np.dot(np.linalg.inv(r), QTb)
    # -->训练集测试结果
    train_result = np.dot(A, x)
    for i in range(len(train_result)):
        if train_result[i] > 0:
            train_result[i] = 1
        else:
            train_result[i] = -1
    train_result.reshape(len(train_result), 1)

    # 5.统计结果
    show(train_labels, train_result, "训练集：")
    # -->测试集测试结果
    A_test = np.hstack([np.ones([10000, 1]), test_images_feature.T])  # 测试集
    b_test = test_labels
    # test_result = A_test.dot(x)
    test_result = np.dot(A_test, x)
    for i in range(len(test_result)):
        if test_result[i] > 0:
            test_result[i] = 1
        else:
            test_result[i] = -1
    test_result.reshape(len(test_result), 1)
    show(test_labels, test_result, "测试集：")
    end_time = time.time()
    print("计算耗时：", end_time - start_time)
